subroutine rgp_fortran_energy_aabb( &
       no, & ! (in)    number of occupied orbitals
       nv, & ! (in)    number of unoccupied orbitals
    Tiajb, & ! (in)    two electron integrals in 1122/aabb order
      Via, & ! (in)    new rows
      Ujb, & ! (in)    new columns
      Pab, & ! (in)    pairing matrix of unoccupied sites (indices are alpha beta)
      Bkl, & ! (in)    inverse matrix (indices are beta alpha)
      Blk, & ! (in)    transpose of inverse matrix (indices are alpha beta)
     VBia, & ! (in)    contraction of new rows with inverse
     BUjb, & ! (in)    contraction of new cols with inverse
    VBUab, & ! (out)   contraction of new rows, new cols, and inverse
     TBab, & ! (out)   contraction of tei and inverse
    TVBjb, & ! (out)   contraction of tei, new rows, and inverse
   retval)   ! (inout) running total for the energy

implicit none

integer,      intent(in)    ::     no
integer,      intent(in)    ::     nv
real(kind=8), intent(in)    ::  Tiajb(no,nv,no,nv)
real(kind=8), intent(in)    ::    Via(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::    Ujb(no,nv) ! second index labels which col
real(kind=8), intent(in)    ::    Pab(nv,nv) ! alpha beta
real(kind=8), intent(in)    ::    Bkl(no,no) ! beta alpha
real(kind=8), intent(in)    ::    Blk(no,no) ! alpha beta
real(kind=8), intent(in)    ::   VBia(no,nv) ! second index labels which row
real(kind=8), intent(in)    ::   BUjb(no,nv) ! second index labels which col
real(kind=8), intent(out)   ::  VBUab(nv,nv)
real(kind=8), intent(out)   ::   TBab(nv,nv)
real(kind=8), intent(out)   ::  TVBjb(no,nv)
real(kind=8), intent(inout) :: retval

integer :: a, b, i, j, k, l

include 'formic/fqmc/rgp_fortran_interface.fpp'

! compute VBUab
VBUab = 0.00d+00
call rgp_fortran_build_vbuab_aabb(no, nv, Via, BUjb, VBUab)

! compute TBab
call rgp_fortran_build_tbab_aabb(no, nv, Tiajb, Blk, TBab)

! compute TVBjb
call rgp_fortran_build_tvbjb(no, nv, Tiajb, VBia, TVBjb)

! term 6 and 16
do b = 1,nv
do a = 1,nv
  retval = retval + TBab(a,b) * ( Pab(a,b) - VBUab(a,b) )
enddo
enddo

! term 12
retval = retval + ddot(no*nv, TVBjb, 1, BUjb, 1)

end subroutine rgp_fortran_energy_aabb
